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Abstract 

We continue our investigation to the use of the variational method to derive 
flow relations for generalized Newtonian fluids in confined geometries. While 
in the previous investigations we used the straight circular tube geometry with 
eight fluid rheological models to demonstrate and establish the variational 
method, the focus here is on the plane long thin slit geometry using those eight 
rheological models, namely: Newtonian, power law, Ree-Eyring, Carreau, 
Cross, Casson, Bingham and Herschel-Bulkley. We demonstrate how the 
variational principle based on minimizing the total stress in the flow conduit 
can be used to derive analytical expressions, which are previously derived by 
other methods, or used in conjunction with numerical procedures to obtain 
numerical solutions which are virtually identical to the solutions obtained 
previously from well established methods of fluid dynamics. In this regard, 
we use the method of Weissenberg-Rabinowitsch-Mooney-Schofield (WRMS), 
with our adaptation from the circular pipe geometry to the long thin slit 
geometry, to derive analytical formulae for the eight types of fluid where these 
derived formulae are used for comparison and validation of the variational 
formulae and numerical solutions. Although some examples may be of little 
value, the optimization principle which the variational method is based upon 
has a significant theoretical value as it reveals the tendency of the flow system 
to assume a configuration that minimizes the total stress. Our proposal also 
offers a new methodology to tackle common problems in fluid dynamics and 
rheology. 

Keywords: Euler-Lagrange variational principle; fluid mechanics; rheology; 
generalized Newtonian fluid; slit flow; pressure-flow rate relation; Newtonian; 
power law; Ree-Eyring; Carreau; Cross; Casson; Bingham; Herschel-Bulkley; 
Weissenberg-Rabinowitsch-Mooney-Schofield method. 
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1 Introduction 


The flow of Newtonian and non-Newtonian fluids in various confined geometries, 
such as tubes and slits, is commonplace in many natural and technological systems. 
Hence, many methods have been proposed and developed to solve the flow problems 
in such geometries applying different physical principles and employing a diverse 
collection of analytical, empirical and numerical techniques. These methods range 
from employing the first principles of fluid dynamics which are based on the rules of 
classical mechanics to more specialized techniques such as the use of Weissenberg- 
Rabinowitsch-Mooney-Schofield relation or one of the Navier-Stokes adaptations 
[ 1 . 2 ]- 

One of the elegant mathematical branches that is regularly employed in the 
physical sciences is the calculus of variation which is based on optimizing functionals 
that describe certain physical phenomena. The variational method is widely used in 
many disciplines of theoretical and applied sciences, such as quantum mechanics and 
statistical physics, as well as many fields of engineering. Apart from its mathematical 
beauty, the method has a big advantage over many other competing methods by 
giving an insight into the investigated phenomena. The method does not only solve 
the problem formally and hence provides a mathematical solution but it also reveals 
the Nature habits and its inclination to economize or lavish on one of the involved 
physical attributes or the other such as time, speed, entropy and energy. Some 
of the well known examples that are based on the variational principle or derived 
from the variational method are the Fermat principle of least time and the curve of 
fastest descent (brachistochrone). These examples, among many other variational 
examples, have played a significant role in the development of the modern natural 
sciences and mathematical methods. 

In reference [3] we made an attempt to exploit the variational method to obtain 
analytic or numeric relations for the flow of generalized Newtonian fluids in confined 
geometries where we postulated that the flow profile in a flow conduit will adjust 
itself to minimize the total stress. This attempt was later [4] extended to include 
more types of non-Newtonian fluids. In the above references, the flow of eight fluid 
models (Newtonian, power law, Bingham, Herschel-Bulkley, Carreau, Cross, Ree- 
Eyring and Casson) in straight cylindrical tubes was investigated analytically and / or 
numerically with some of these models confirming the stated variational hypothesis 
while others, due to mathematical difficulties or limitation of the underlying principle, 
demonstrated behavioral trends that are consistent with the variational hypothesis. 

No mathematically rigorous proof was gives in [3, 4] to establish the proposed 
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variational method that is based on minimizing the total stress in its generality. 
Furthermore, we do not make any attempt here to present such a proof. However, 
in the present paper we make an attempt to consolidate our previous proposal and 
findings by giving more examples, this time from the slit geometry rather than 
the tube geometry, to validate the use of the variational principle in deriving flow 
relations in confined geometries for generalized Newtonian fluids. 

The plan for this paper is as follow. In the next section § 2 we present the 
general formulation of the variational method as applied to the long thin slits and 
derive the main variational equation that will be used in obtaining the flow relations 
for the generalized Newtonian fluids. In section § 3 we apply and validate the 
variational method for five types of non-viscoplastic fluids, namely: Newtonian, 
power law, Ree-Eyring, Carreau and Cross; while in section § 4 we apply and validate 
the method for three types of viscoplastic fluids, namely: Casson, Bingham and 
Herschel-Bulkley. We separate the viscoplastic from the non-viscoplastic because 
the variational method strictly applies only to non- viscoplastic fluids, and hence its 
use with viscoplastic fluids is an approximation which is valid and good only when 
the yield stress value of these fluids is low and hence the departure from fluidity is 
minor. In the validation of both non-viscoplastic and viscoplastic types we use the 
aforementioned WRMS method where we compare the variational solutions to the 
analytical solutions obtained from the WRMS method where analytical formulae 
for the eight types of fluid are derived in the Appendix (§ 7). The paper in ended in 
section § 5 where general discussion and conclusions about the paper, its objectives 
and achievements are presented. 

2 Method 

The rheological behavior of generalized Newtonian fluids in one dimensional shear 
flow is described by the following constitutive relation 


t = F7 (!) 

where r is the shear stress, 7 is the rate of shear strain, and /i is the shear viscosity 
which is normally a function of the contemporary rate of shear strain but not of the 
deformation history although it may also be a function of other physical parameters 
such as temperature and pressure. The latter parameters are not considered in 
the present investigation as we assume a static physical setting (i.e. isothermal, 
isobaric, etc.) apart from the purely kinematical aspects of the deformation process 
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that is necessary to initiate and sustain the flow. 

In the following we use the slit geometry, depicted in Figure 1, as our flow 
apparatus where 2 B is the slit thickness, L is the length of the slit across which 
a pressure drop A p is imposed, and W is the part of the slit width that is under 
consideration although for the purpose of eliminating lateral edge effects we assume 
that the total width of the slit is much larger than the considered part W. We 
choose our coordinates system so that the slit smallest dimension is being positioned 
symmetrically with respect to the plane z = 0. 



Figure 1: Schematic drawing of the slit geometry which is used in the present 
investigation. 

For the slit geometry of Figure 1 the total stress is given by 


T t = 


f-r+s 


'T-B 


dr = r d ^dz = 
J -B dz 


+B d 

dz (n)dZ = 


" +B 'dll d'y' 

7 -jz +^)dz (2) 


l-B 


l-B 


dz dz 


where Tt is the total stress, and t±b is the shear stress at the slit walls corresponding 
to z = ±B. 

The total stress, as given by Equation 2, can be minimized by applying the 
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Euler-Lagrange variational principle which, in one of its forms, is given by 


A f f _ d L = n 

dx \ ^ dy' ) dx 

where the symbols corresponding to our problem statement are defined as 


( 3 ) 


dy d^ df_dfdy dy/\ 

I = 2 ’ V = J ’ / = hfe +fl d? alld a^ = ay {^. + ^) =11 (4) 

On substituting these symbols into Equation 3 the following equation is obtained 


d f dy d'y dj\ d f dy dy' 
iz yTz + ^Tz ~ ^Tz ) ~ Tz yTz + 1 = 0 

Considering the fact that for the considered flow systems 

dy d'y 

7 — + /r— = G 
dz dz 


(5) 


( 6 ) 


where G is a constant, it can be shown that Equation 5 can be reduced to two 
independent variational forms 


d_ ( dy 

dz \ dz * 


(7) 


and 


d ( d l\ _ n 

dz Vdz * 


( 8 ) 


In the following two sections we use the second form where we outline the 
application of the variational method, as summarized in Equation 8 , to validate and 
demonstrate the use of the variational method to obtain flow relations correlating 
the volumetric flow rate Q through the slit to the pressure drop A p across the 
slit length L for generalized Newtonian fluids. The plan is that we derive fully 
analytical expressions when this is viable, as in the case of Newtonian and power 
law fluids, and partly analytical solutions when the former is not viable, e.g. for 
Ree-Eyring and Herschel-Bulkley fluids. In the latter case, the solution is obtained 
numerically in its final stages, following a variationally based derivation, by using 
numerical integration and simple numerical solvers. 

In this investigation we assume a laminar, incompressible, time-independent, 
fully-developed, isothermal flow where entry and exit edge effects are negligible. We 
also assume negligible body forces and a blunt flow speed profile with a no-shear 
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stationary region at the profile center plane which is consistent with the considered 
type of fluids and flow conditions, i.e. viscous generalized Newtonian fluids in a 
pressure-driven laminar flow. As for the plane slit geometry, we assume, following 
what is stated in the literature, a long thin slit with B <C IT and B <C L although we 
believe that some of these conditions are redundant according to our own statement 
and problem settings. We also assume that the slit is rigid and uniform in shape 
and size, that is its walls are not made of deformable materials, such as elastic or 
viscoelastic, and the slit does not experience an abrupt or gradual change in B. 

3 Non-Viscoplastic Fluids 

For non- viscoplastic fluids, the variational principle strictly applies. In this section 
we apply the variational method to five non-viscoplastic fluids and compare the 
variational solutions to the analytical solutions obtained from the WRMS method. 
These fluids are: Newtonian, power law, Ree-Eyring, Carreau and Cross. All these 
models have analytical solutions that can be obtained from various traditional 
methods of fluid dynamics which are not based on the variational principle. Hence 
the agreement between the solutions obtained from the traditional methods with 
the solutions obtained from the variational method will validate and vindicate the 
variational approach. 

As indicated earlier, to derive non-variational analytical relations we use a 
method similar to the one ascribed to Weissenberg, Rabinowitsch, Mooney, and 
Schofield [1] for the flow of generalized Newtonian fluids in uniform tubes with 
circular cross sections, where we adapt and apply the procedure to long thin slits, and 
hence we label this method with WRMS to abbreviate the names of its originators. 
The WRMS method is fully explained and applied in the Appendix (§ 7) to derive 
analytical relations to all the eight fluid models that are used in the present paper. 
For some of these fluids full analytical solutions from the variational principle 
are obtained and hence a direct comparison between the analytical expressions 
obtained from the two methods can be made, while for other fluids a mixed 
analytical and numerical procedure is employed to obtain numerical solutions from 
the variational principle and hence a representative sample of numerical solutions 
from both methods is presented for comparison and validation, as will be clarified 
and demonstrated in the following subsections. 
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3.1 Newtonian 

The viscosity of Newtonian fluids is constant, that is 


d do 

and therefore Equation 8 becomes 


( 9 ) 


d 

dz 



0 


On performing the outer integration we obtain 


( 10 ) 


d°2~ A ( n ) 

where A is the constant of integration. On performing the inner integration we 
obtain 


7 = —z + D (12) 

do 

where id is a second constant of integration. Now from the no-shear condition at 
the slit center plane z = 0, D can be determined, that is 


7 (~ = 0) = 0 => D = 0 (13) 

Similarly, from the no-slip boundary condition [5] at z — AB which controls the 
wall shear stress we determine A, i.e. 

F± 2BW Ap BAp 

T±B = ^ = = -r (14) 

where t±b is the shear stress at the slit walls, Fj_ is the flow driving force which is 
normal to the slit cross section in the flow direction, and cr|| is the area of the slit 
walls which is tangential to the flow direction. Hence 


Therefore 


7 (z = ±B) 


t±b 

do 


BAp AB 

doL do 


A p 

1 ~T 


(15) 


7(z) 


A p 

— T* 
fJ'o-LJ 


(16) 


On integrating the rate of shear strain with respect to z, the standard parabolic 
speed profile is obtained, that is 
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where v(z) is the fluid speed at 0 in the x direction and E is another constant of 
integration which can be determined from the no-slip at the wall boundary condition, 
that is 


i.e. 


v (z = ±B ) = 0 


=> E 


-PLb* 


(18) 




(19) 


The volumetric flow rate is then obtained by integrating the flow speed profile over 
the slit cross sectional area in the z direction, that is 


f +b 


Q = 


vWdz 


l-B 


2W Ap 
2/r 0 T 



WAp 

p 0 L 




( 20 ) 


that is 


0 = ( 21 ) 

which is the well known volumetric flow rate formula for the flow of Newtonian 
fluids in a plane long thin slit as obtained by other methods which are not based on 
the variational principle. This formula is derived in the Appendix (§ 7, Equation 
80) using the WRMS method. It also can be found in several classic textbooks of 
fluid mechanics, e.g. Bird et al. [2| Table 4.5-14 where p a = p and Ap = P 0 — P L . 


3.2 Power Law 

The shear dependent viscosity of power law fluids is given by [1, 2, 6] 

(Jt = ki 1 - 1 (22) 

where k is the power law viscosity consistency coefficient and n is the flow behavior 
index. On applying the Euler-Lagrange variational principle (Equation 8) we obtain 
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4 (*■"£)- 

On performing the outer integral we obtain 

k 1 n ~ 1 — = A (24) 

dz 

On separating the two variables in the last equation and integrating both sides we 
obtain 


7 = ^ ( Az + D ) (25) 

where A and D are the constants of integration which can be determined from the 
two limiting conditions, that is 


7 (z = 0) = 0 


D = 0 


(26) 


and 


7 (z = B)= {/y = 


BA P _ JVlar 
L k V k A 


< 27 > 


where the first step in the last equation is obtained from the constitutive relation 
of power law fluids, i.e. 


r = k^ n (28) 

with the substitution z = B in Equation 25. Hence, from Equation 25 we obtain 

< 29 > 

On integrating the rate of shear strain with respect to z, the flow speed profile is 
determined, i.e. 


v(z) 



= 

kL 


n 

n + 1 


^z 1+1/n +E 

kL 


(30) 

where E is another constant of integration which can be determined from the no-slip 
at the wall condition, that is 
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v (z — B) — 0 =>■ E 


(31) 


= _0_ ,”/A> Rl+l/r. 

n + 1 V fcL 

i.e. 


"W = ^i\l( Bl+,/ "- 1+I/ “) (32) 

The volumetric flow rate can then be obtained by integrating the flow speed profile 
with respect to the cross sectional area in the z direction, that is 


i.e. 


f +B 2 Wn 

Q= v Wdz = - 

J-B n + 1 




2 Wn n l Ap 

n + 1 V kL 




z 2+l/n 

2 + l/n_ 




2iTn n lAp 

n + 1 V fcL 


j^2+l/n 


j^2+l/n ' 

2 + l/n 


(33) 

(34) 

(35) 


2HT;" j h n /BAp 

- _ 2n + l V feL 


(36) 


which is the well known volumetric flow rate relation for the flow of power law fluids 
in a long thin slit as obtained by other non-variational methods. This formula is 
derived in the Appendix (§ 7, Equation 84) using the WR.MS method. It can also 
be found in textbooks of fluid mechanics such as Bird et al. [2| Table 4.2-1 where 
k = m and Ap = P 0 — Pi. 


3.3 Ree-Eyring 

For Ree-Eyring fluids, the constitutive relation between shear stress and rate of 
strain is given by [2] 

r = r c asinh (37) 

where r c is a characteristic shear stress and /i r is the viscosity at vanishing rate of 
strain. Hence, the generalized Newtonian viscosity is given by 
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T r c asinh 



( 38 ) 


7 


7 


On substituting /i from the last relation into Equation 8 we obtain 



dz 1 7 


0 


(39) 


On integrating once we get 


t c asinh 



A 


(40) 


7 


dz 


where A is a constant. On separating the two variables and integrating again we 
obtain 


where the constant of integration D is absorbed within the integral on the left hand 
side. The integral on the left hand side of Equation 41 when evaluated analytically 
produces an expression that involves logarithmic and polylogaritlnnic functions 
which when computed produce significant errors especially in the neighborhood of 
z — 0. To solve this problem we used a numerical integration procedure to evaluate 
this integral, and hence obtain A, numerically using the boundary condition at the 
slit wall, that is 


where Tb is given by Equation 14. This was then followed by obtaining 7 as a 
function of z using a bisection numerical solver in conjunction with a numerical 
integration procedure based on Equation 41. Due to the fact that the constant 
of integration, D , is absorbed in the left hand side and a numerical integration 
procedure was used rather than an analytical evaluation of the integral on the left 
hand side of Equation 41, there was no need for an analytical evaluation of this 
constant using the boundary condition at the slit center plane, i.e. 


/ 



(41) 


7 



(42) 


7 (z = 0 ) = 0 


(43) 
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The numerically obtained 7 was then integrated numerically with respect to z to 
obtain the flow speed as a function of z where the 110 -slip boundary condition at 
the wall is used to have an initial value v = 0 that is incremented on moving inward 
from the wall toward the center plane. The flow speed profile was then integrated 
numerically with respect to the cross sectional area to obtain the volumetric flow 
rate. 

To test the validity of the variational method we made extensive comparisons 
between the solutions obtained from the variational method to those obtained from 
the WRMS method using widely varying ranges of fluid and slit parameters. In 
Figure 2 we compare the WRMS analytical solutions as derived in the Appendix (§ 
7, Equation 88 ) with the variational solutions for two sample cases. As seen, the two 
methods agree very well which is typical in all the investigated cases. The minor 
differences between the solutions of the two methods can be easily explained by the 
accumulated errors arising from repetitive use of numerical integration and numerical 
solvers in the variational method. The errors as estimated by the percentage relative 
difference are typically less than 0.5% when using reliable numerical integration 
schemes with reasonably fine discretization mesh and tiny error margin for the 
convergence condition of the numerical solver. This is also true in general for the 
other types of fluid that will follow in this section. 

3.4 Carreau 

For Carreau fluids, the viscosity is given by [2, 6 - 8 ] 


where po is the zero-shear viscosity, /i,- is the infinite-shear viscosity, A is a characteris- 
tic time constant, and n is the flow behavior index. On applying the Euler-Lagrange 
variational principle (Equation 8 ) and following the derivation, as outlined in the 
previous subsections, we obtain 


where 2 Ti is the hypergeometric function of the given argument with its real part 
being used in the computation, and A and D are the constants of integration. From 
the two boundary conditions at z = 0 and z = B, A and D can be determined, that 

is 


h ~ hi + (ho — hi) [l + -^ 2 7 2 ] 


2Jl("— 1)/2 


(44) 



(45) 
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(a) (b) 


Figure 2: Comparing the WRMS analytical solutions, as given by Equation 88, to 
the variational solutions for the flow of Ree-Eyring fluids in long thin slits with (a) 
fi r = 0.005 Pa.s, t c = 600 Pa, B = 0.01 m, W = 1.0 m and L = 1.0 m; and (b) 
p r = 0.015 Pa.s, r c = 400 Pa, B = 0.013 m, W = 1.0 nr and L = 1.7 m. In both 
sub-figures, the vertical axis represents the volumetric flow rate, Q, in m 3 .s -1 while 
the horizontal axis represents the pressure drop, A p. in Pa. The average percentage 
relative difference between the WRMS analytical solutions and the variational 
solutions for these cases are about 0.38% and 0.39% respectively. 


7 (~ = 0) = 0 =*► D = 0 (46) 

and 


7 (z = B) 
where y# is 


= 7 b 

the shear rate 


, , „ , 1 1 — n 

VilB + (hO — hi) 7 B 2-f 1 ( 2 ’ 


at the slit wall. Now, by definition 



we have 


AB 

(47) 


hn 7 s — t b 


(48) 


that is 


hi + (ho - hi) [l + A 2 y|] (n 1)/2 


7 b 


BAp 

L 


(49) 


From the last equation, y# can be obtained numerically by a numerical solver, based 
for example on a bisection method, and hence from Equation 47 A is obtained. 
Equation 45 can then be solved numerically to obtain the shear rate y as a function 
of z. This will be followed by integrating y numerically with respect to z to obtain 
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the speed profile, v(r), where the no-slip boundary condition at the wall is used to 
have an initial value v(z = ±R) = 0 that is incremented on moving inward toward 
the center plane. The speed profile, in its turn, will be integrated numerically with 
respect to the slit cross sectional area to obtain the volumetric flow rate Q. 

In Figure 3 we present two sample cases for the flow of Carreau fluids in thin 
slits where the WRMS analytical solutions, as given by Equation 96, are compared 
to the variational solutions. Good agreement can be seen in these plots which 
are typical of the investigated cases. The main reason for the difference between 
the WRMS and variational solutions is the persistent use of numerical solvers and 
numerical integration in the implementation of the variational method as well as the 
use of the hypergeometric function in both methods. The numerical implementation 
of this function can cause instability and failure to converge satisfactorily in some 


cases. 


3.5 Cross 

For Cross fluids, the viscosity is given by [ 6 , 9] 


H — Ah + 


HO Hi 


(50) 


1 + \m^m 

where ft 0 is the zero-shear viscosity, /i* is the infinite-shear viscosity, A is a character- 
istic time constant, and m is an indicia! parameter. Following a similar derivation 
method to that outlined in Carreau, we obtain 


I 777 “I - 1 

Hi 7 + (ho - Hi) 7 2 C ( 1, — ; ; -A m 7 m ) = Az 

m m 


(51) 


where 


A = 


M<7b + (Mo - 7 b 2 Fi (1, Jj; =£1; -A”> 7 £) 


B 


(52) 


with 7 b being obtained numerically as outlined in Carreau. Equation 51 can then 
be solved numerically to obtain the shear rate 7 as a function of z. This is followed 
by obtaining v from 7 and Q from v by using numerical integration, as before. 

In Figure 4 we present two sample cases for the flow of Cross fluids in thin slits 
where we compare the WRMS analytical solutions, as given by Equation 104, to the 
variational solutions. As seen in these plots, the agreement is good with the main 
reason for the departure between the two methods is the persistent use of numerical 
solvers and numerical integration in the variational method as well as the use of 
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Figure 3: Comparing the WRMS analytical solutions, as given by Equation 96, to 
the variational solutions for the flow of Carreau fluids in long thin slits with (a) 
n = 0.9, fi 0 = 0.13 Pa.s, /i* = 0.002 Pa.s, A = 0.8 s, B = 0.011 m, W = 1.0 m 
and L = 1.25 m; and (b) n = 0.85, /io = 0.15 Pa.s, /ij = 0.01 Pa.s, A = 1.65 s, 
B = 0.011 m, W = 1.0 m and L — 1.4 m. The axes are as in Figure 2, while the 
differences are about 0.21% and 0.25% respectively. 


the hypergeometric function in both methods, as explained in the case of Carreau. 


4 Viscoplastic Fluids 

The yield stress fluids are not strictly subject to the variational principle due to 
the existence of a solid non-yield region at the center which does not obey the 
stress optimization condition and hence the Euler-Lagrange variational method 
is not strictly applicable to these fluids. However, the method provides a good 
approximation when the value of the yield stress is low so that the effect of the 
11011 -yield region at and around the center plane of the slit on the flow profile is 
minor. In the following subsections we apply the variational method to three yield 
stress fluids and obtain some solutions from sample cases which are representative 
of the many cases that were investigated. 


4.1 Casson 

For Casson fluids, the constitutive relation is given by [2, 6] 


r = 


(A'7) 1/2 + t 0 1/2 1 


(53) 


where K is the viscosity consistency coefficient, and r D is the yield stress. Hence 
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Figure 4: Comparing the WRMS analytical solutions, as given by Equation 104, 
to the variational solutions for the flow of Cross fluids in long thin slits with (a) 
m = 0.65, /ro = 0.032 Pa.s, /q = 0.004 Pa.s, A = 4.5 s, B = 0.013 m, W = 1.0 m 
and L = 1.3 m; and (b) m = 0.5, /io = 0.075 Pa.s, /q = 0.004 Pa.s, A = 0.8 s, 
B = 0.006 m, W = 1.0 m and L = 0.8 m. The axes are as in Figure 2, while the 
differences are about 0.29% and 0.88% respectively. 


r 

h = - = 

7 


(A' 7 ) 1/2 + tP 


7 


(54) 


On substituting /i from the last equation into the main variational relation, as given 
by Equation 8, we obtain 


d 

dz 


7 dz 

On integrating twice with respect to z we obtain 


(Kl) 1 ' 2 + r'J 2 


2 ^ 


= 0 


(55) 


K'y + 4 (Kt q 7) 1//2 + t 0 In (7) = Az + D 


(56) 


where A and D are constants. Now, from the boundary condition at the slit center 
plane we have 


7 (z = 0) = 0 (57) 

so we set D = 0 to constrain the solution at z — 0. As for the second boundary 
condition at the slit wall, z = B, we have 
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K 7 b + 4 (Kt 0 ib) 1/2 + t 0 In ( 75 ) = AB 


(58) 


where the rate of shear strain at the slit wall, 7 #, is obtained from applying the 
constitutive relation at the wall and hence is given by 


7 b = 


A /TB To 


1/2 


Ji 


(59) 


with the shear stress at the slit wall, 7 #, being given by Equation 14. Hence 


_ K 1b + 4 (Jl t 0 7 b ) 1/2 + t 0 In ( 75 ) 

B 

Equation 56 defines 7 implicitly in terms of 0 and hence it is solved numerically 
using, for instance, a numerical bisection method to find 7 as a function of z with 
avoidance of the very immediate neighborhood of z = 0 which, as explained earlier, 
is not subject to the variational method. This is equivalent to integrating between 
r 0 and t b , rather than between 0 and t b , in the WRMS method, as employed in 
the Appendix (§ 7), for the case of Casson, Bingham and Herschel-Bulkley fluids. 
Although the value of z that defines the start of the yield region near the center is 
not known exactly, we already assumed that the use of the variational method is only 
legitimate when r 0 is small and hence the non- yield plug region is small and hence 
its effect is minor, so any error from an ambiguity in the exact limit of the integral 
near the z = 0 should be negligible especially at high flow rates where this region 
shrinks and hence using a lower limit of the integral at the immediate neighborhood 
of z = 0 will give a more exact definition of the yield region. On obtaining 7 
numerically, v and 0 can be obtained successively by numerical integration, as 
before. 

In Figure 5 we present two sample cases for the flow of Casson fluids in thin slits 
where we compare the WRMS analytical solutions, as given by Equation 110, with 
the variational solutions. As seen, the agreement is reasonably good considering 
that the variational method is just an approximation and hence it is not supposed 
to produce identical results to the analytically derived solutions. The two plots 
also indicate that the approximation is worsened as the value of the yield stress 
increases, resulting in the increase of the effect of the non-yield region at the center 
plane of the slit which is not subject to the variational principle, and hence more 
deviation between the two methods is observed. 
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Figure 5: Comparing the WRMS analytical solutions, as given by Equation 110, 
to the variational solutions for the flow of Casson fluids in long thin slits with (a) 
K = 0.025 Pa.s, r D = 0.1 Pa, B = 0.01 m, W = 1.0 m and L = 0.5 m; and (b) 
K = 0.05 Pa.s, t 0 = 0.5 Pa, B = 0.025 m, W = 1.0 m and L = 1.5 m. The axes are 
as in Figure 2, while the differences are about 0.57% and 2.84% respectively. 

4.2 Bingham 

For Bingham fluids, the viscosity is given by [1, 2, 6] 

ii = - + a (6i) 

7 

where r a is the yield stress and C is the viscosity consistency coefficient. On 
applying the variational principle, as formulated by Equation 8, and following the 
previously outlined method we obtain 


t 0 In 7 + C ,r j = Az + D (62) 

where A and D are the constants of integration. Using the boundary conditions at 
the center plane and at the slit wall and following a similar procedure to that of 
Casson, we obtain 


D = 0 


and 


A 


r ° i 

B 11 





(63) 


The strain rate is then obtained numerically from Equation 62, and thereby v and 
Q are computed successively, as explained before. 

In Figure 6 two sample cases of the flow of Bingham fluids in thin slits are 
presented. As seen, the agreement between the WRMS solutions, as obtained from 
Equation 113, and the variational solutions are rather good despite the fact that 
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the variational method is an approximation when applied to viscoplastic fluids. 


4.3 Herschel-Bulkley 

The viscosity of Herschel-Bulkley fluids is given by [1, 2, 6 ] 


n = ~ + cy 1 - 1 (64) 

7 

where r D is the yield stress, C is the viscosity consistency coefficient and n is the 
flow behavior index. On following a procedure similar to the procedure of Bingham 
model with the application of the 7 two boundary conditions, we get the following 
equation 


where 


C 

r 0 In 7 -| 7 ” = Az 

n 


(65) 


A 


— In 

B \LC 






( 66 ) 


On solving Equation 65 numerically, 7 as a function of z is obtained, followed by 
obtaining v and 0. as explained previously. 

In Figure 7 we compare the WRMS analytical solutions of Equation 117 to the 
variational solutions for two sample Herschel-Bulkley fluids, one shear thinning 
and one shear thickening, both with yield stress. As seen, the agreement between 
the solutions of the two methods is good as in the previous cases of Casson and 
Bingham fluids. 


5 Conclusions 

In this paper we provided further evidence for the validity of the variational 
method which is based on minimizing the total stress in the flow conduit to 
obtain flow relations for the generalized Newtonian fluids in confined geometries. 
Our investigation in the present paper, which is related to the plane long thin 
slit geometry, confirms our previous findings which were established using the 
straight circular uniform tube geometry. Eight fluid types are used in the present 
investigation: Newtonian, power law, Ree-Eyring, Carreau, Cross, Casson, Bingham 
and Herschel-Bulkley. This effort, added to the previous investigations, should be 
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Figure 6: Comparing the WRMS analytical solutions, as given by Equation 113, to 
the variational solutions for the flow of Bingham fluids in long thin slits with (a) 
C = 0.02 Pa.s, t 0 = 0.25 Pa, B = 0.015 nr, W = 1.0 m and L = 0.75 m; and (b) 
C = 0.034 Pa.s, t 0 = 0.75 Pa, B = 0.018 m, W = 1.0 m and L = 1.25 m. The axes 
are as in Figure 2, while the differences are about 1.12% and 1.96% respectively. 




Figure 7: Comparing the WRMS analytical solutions, as given by Equation 117, to 
the variational solutions for the flow of Herschel-Bulkley fluids in long thin slits with 
(a) n = 0.8, C = 0.05 Pa.s n , r G = 0.5 Pa, B = 0.01 m, W = 1.0 m and L = 1.2 m; 
and (b) n = 1.25, C = 0.025 Pa.s n , r 0 = 1.0 Pa, B = 0.03 m, W = 1.0 m and 
L = 1.3 m. The axes are as in Figure 2, while the differences are about 1.49% and 
1.74% respectively. 
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sufficient to establish the variational method and the optimization principle upon 
which the method rests. For the Newtonian and power law fluids, full analytical 
solutions are obtained from the variational method, while for the other fluids mixed 
analytical-numerical procedures were established and used to obtain the solutions. 

Although some of the derived expressions and solutions are not of interest of 
their own as they can be easily obtained from other non-variational methods, the 
theoretical aspect of our investigation should be of great interest as it reveals 
a tendency of the flow system to minimize the total stress which the variational 
method is based upon: hence giving an insight into the underlying physical principles 
that control the flow of fluids. 

The value of our investigation is not limited to the above mentioned theoretical 
aspect but it has a practical aspect as well since the variational method can be 
used as an alternative to other methods for other geometries and other rheological 
fluid models where mathematical difficulties may be overcome in one formulation 
based on one of these methods but not the others. The variational method is also 
more general and hence enjoys a wider applicability than some of the other methods 
which are based on more special or restrictive physical or mathematical principles. 
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6 


Nomenclature 


7 rate of shear strain (s 

5 Ho- ^ (Pa.s) 

A characteristic time constant (s) 

H fluid shear viscosity (Pa.s) 

Ho zero-shear viscosity (Pa.s) 

Hi infinite-shear viscosity (Pa.s) 

Ho Newtonian viscosity (Pa.s) 

H r low-shear viscosity in Ree-Eyring model (Pa.s) 

<7 1| area of slit wall tangential to the flow direction (m 2 ) 
r shear stress (Pa) 

t±b shear stress at slit walls corresponding to z = ±P (Pa) 
r c characteristic shear stress in Ree-Eyring model (Pa) 
t 0 yield stress in Casson, Bingham and Herschel-Bulkley models (Pa) 
r t total shear stress (Pa) 

B slit half thickness (m) 

C viscosity consistency coefficient in Herschel-Bulkley model (Pa.s n ) 

C" viscosity consistency coefficient in Bingham model (Pa.s) 

/ A m 7 £ 

2 -Fi hypergeometric function 

F± force normal to the slit cross section (N) 

9 1 + / 

Jca definite integral expression for Carreau model (Pa 2 .s _1 ) 

Jcr definite integral expression for Cross model (Pa 2 . s' 1 ) 
k viscosity consistency coefficient in power law model (Pa.s n ) 

K viscosity consistency coefficient in Casson model (Pa.s) 

L slit length (m) 

m indicia! parameter in Cross model 

n flow behavior index in power law, Carreau and Herschel-Bulkley models 
A p pressure drop across the slit length (Pa) 

Po pressure at the slit entrance (Pa) 

Pl pressure at the slit exit (Pa) 
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Q volumetric flow rate (m 3 .s -1 ) 
v fluid speed in the flow direction (m.s -1 ) 
W slit width (m) 

z coordinate of slit smallest dimension (m) 
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7 Appendix 

Here we derive a general formula for the volumetric flow rate of generalized Newto- 
nian fluids in rigid plane long thin uniform slits following a method similar to the 
one ascribed to Weissenberg, Rabinowitsch, Mooney, and Schofield [1], and hence 
we label the method with WRMS. We then apply it to obtain analytical relations 
correlating the volumetric flow rate to the pressure drop for the flow of Newtonian 
and seven non-Newtonian fluids through the above described slit geometry. 

The differential flow rate through a differential strip along the slit width is given 

by 


dQ = vWdz 


(67) 


where Q is the volumetric flow rate and v = v(z) is the fluid speed at z in the x 
direction according to the coordinates system demonstrated in Figure 1. Hence 


Q_ 

w 


f +B 


vdz 


( 68 ) 


f-B 


On integrating by parts we get 


Q r i+S 

w = 


w+B 


zdv 


(69) 


'V_ B 


The first term on the right hand side is zero due to the no-slip boundary condition, 
and hence we have 


Q_ 

W 


f V+B 


zdv 


IV_ B 


Now, from Equation 14, we have 


t±b = 


BAp 


where t±b is the shear stress at the slit walls. Similarly we have 

zAp 


r z = 


where t 2 is the shear stress at z. Hence 


z Bt z 

T z = => Z= 

B t± b 

We also have 


dz = 


Bdr z 

t±b 


(70) 


(71) 


(72) 


(73) 
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lz = 


dv 

dz 


, , Bdr z 

dv = -7 Z dz = —'-fz 

t±b 


(74) 


Now due to the symmetry with respect to the plane z = 0 we have 


t b = r +B = t_ b (75) 

On substituting from Equations 73 and 74 into Equation 70, considering the flow 
symmetry across the center plane z = 0, and changing the limits of integration we 
obtain 


that is 


Q_ 

W 


rr+B 


' T —B 


Bt z 

t±b 


-lz 


Bdr z 

r±B 



(76) 



where it is understood that 7 = = 7 (z) and r = r z = t(z). 


For Newtonian fluids with viscosity fj, a we have 


r 


T = Hoi 


1 = 


Hence 


( B 


rrs 


1^0 


2W ( b \ 2 r TB 


Q = 2 W I — ) / 7 rdr = ‘ Z1 — ( — ) / r 2 dr = 

\ t bJ Jo Ho \ t bJ Jo 2 !/i d 


2WB" 


that is 


(77) 


(78) 


tb (79) 


Q = 


2WB s Ap 

3/Jj 0 L 


(80) 


For power law fluids we have 

r = /c7 n 


'T\ l / n 

^ 7 = [k) 


Hence 


(81) 


B\ 2 f TB (r\ V' 


tb) 


2W ( B\ 


Q = 2W I — 1 / (-) rdr = ( — ) / r 1+1/n dr 


kJ 


k 1 / n \tb J 


"T B 


(82) 
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Q = 


2 W 


k 1 / 11 (2 + 1/n) 


t bJ 


T. 


2+1/ra 


2WB' 1 


_l/n 


B 


ki/n (2 + 1/n) 


'B 


(83) 


that is 


^ 2 W&n „[BE p 

| q = dtttViyI (84) 

When n — 1, with k = fi 0 , the formula reduces to the Newtonian, Equation 80, as 
it should be. 


For Ree-Eyring fluids we have 


Hence 


that is 


t = T r asinh 


Hrl 


=> 7 = — sinh ( - 

Hr \T C 


Q = 


2W r r ( B\ 


Hr V r B 7 


- 


r 


rsinh — dr 


T r 


Q = 


2 Wr 2 c ( B\ 


2 r 


Hr \t b ) 


t cosh ( — ) — r c sinh ( — 
T r 1 \ T r 


tb 


(85) 


( 86 ) 

(87) 


Q ^2Wrl(By 

Hr \ TB J 


tb cosh ( — j — t c sinh ( — 


( 88 ) 


For Carreau fluids, the viscosity is given by 


H — ~ — Ah + $ (l + A 2 7 2 ) 

7 

where 5 = (jM) — Hi) and n! = (n — 1). Therefore 


2 2\ ri ' /2 


(89) 


r = 7 


2 „, 2\ n ’/ 2 


Hi + $ (l + A 2 7 2 ) 


and hence 


(90) 


dr = 


A, + a (i + AY)" 72 + h'aaY (i + aV ) 1 "- 2,/2 <i 7 (9i) 


2 _, 2 \ O '— 2)/2 


Now, from the WRMS method we have 
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wh = SJ 1TdT (92) 

If we label the integral on the right hand side of Equation 92 with Jc a and substitute 
for r from Equation 90, substitute for dr from Equation 91, and change the 
integration limits we obtain 


riB 


I Ca — 


7 


m + s (i + aV)"' /2 1 U + s (i + a 2 t 2 )” 72 + n'av (i + a 2 7 2 ) ( ”'- 2)/2 


(93) 

On solving this integral equation analytically and evaluating it at its two limits we 
obtain 


dq 


I Ca — 


+ 

+ 


n' 5 2 7 B [ 2 Fi (§, 1 - n'\ §; -A 2 q|) - 2 ^i (§, -n'; |; -A 2 q|)] 


A 2 


(1 + m!) <S 2 q| 2 Fi (|, — n 7 ; |; -A 2 q|) 

3 

n'^i7B [2^1 (|, 1 - y; |; -a 2 7 2 b) - 2 f 1 (|, -f , !; -a 2 7 |)] 

A 2 

(2 + n') 5/Xi7| 2-Pi (|, - y; |; -A 2 7l) + /^7s 


(94) 


where 2 Ti is the hypergeometric function of the given arguments with its real part 
being used in this evaluation. Now, from applying the rheological equation at the 
slit wall we have 


22 \ n '/ 2 


fii + S (l + A 2 7b) 


7 b = 


B A p 


(95) 


From the last equation, 7 b can be obtained numerically by a simple numerical 
solver, such as bisection, and hence Jc a is computed. The volumetric flow rate is 
then obtained from 


Q 


2WB 2 I Ca 


For Cross fluids, the viscosity is given by 


h 


r 

7 


5 

+ 1 + A m 7 m 


(96) 


(97) 
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where 5 = (/i 0 — Hi)- Therefore 


and hence 


r = 7 



* ) 
1 + \m^mj 


(98) 


, / 5 rn<5A m 7 m \ , 

dr = /ij + 5 - d'y (99) 

V 1 + A m 7 m (l + A m 7 m )V 

If we follow a similar procedure to that of Carreau by applying the WRMS method 
and labeling the right hand side integral with A>, substituting for r and dr from 
Equations 98 and 99 respectively, and changing the integration limits of /<> we get 


07 B 


ICr ~ 


7 + 


1 + A m 7 r 


Hi + 


5 


m5\ m Y‘ 


1 + A m 7 m (1 + \mrymy 


d'y ( 100 ) 


On solving this integral equation analytically and evaluating it at its two limits we 
obtain 


= 


[35 2 (m — g) — {<5 2 (m - 3) + 2 mSfu} g 2 2 Fi (l, 1 + -/) + 


Qmg 2 


( 101 ) 


where 


/ = A m 7 B , 9 = 1 + f (102) 

and 2 Ei is the hypergeometric function of the given argument with its real part 
being used in this evaluation. As before, from applying the rheological equation at 
the wall we have 


Hi + 


5 


7 b — 


B A p 


(103) 


1 + A m 7^ 

From this equation, 7 b can be obtained numerically and hence /<> is computed. 
Finally, the volumetric flow rate is obtained from 


Q 


2WB 2 I Cr 


( 104 ) 


For Casson fluids we have 
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t 1/2 = (A'q) 1/2 + r] /2 (r > r 0 ) 


Hence 


7 = 


r^_ r y^ 2 


A' 


On applying the WRMS equation we get 




2 W f BY 


rr B 


Q = [^rj J ( r2 _ c Z\[T r oT i/2 + r 0 r) dr 


that is 


e-¥® 


7 y _ 4 v /7yr 5 / 2 r 0 r 2 

3 5 2 






_ 4/^4 /2 vr| _ 7^ 
3 5 2 30 


(105) 

(106) 

(107) 

(108) 

(109) 

( 110 ) 


For Bingham fluids we have 

T = To + C'f => 7 = ~qJ~ ( r - T °) (m) 

Hence 


Q 


2W 

~cF 



r 0 ) dr 


2 W 

~cy~ 



1~B 


To 


that is 


( 112 ) 


Q 


2W(B\ 2 \t% t q t 2 rf 

O Vr B J [3 2 6 _ 


(113) 


When r Q = 0, with C = /i Q , the formula reduces to the Newtonian, Equation 80, as 
it should be. 
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For Herschel-Bulkley fluids we have 


t = r a + C 7 n 




(r > r 0 ) (114) 


Hence 


that is 


Q 


2 W f B\ 
C^Vb) 




Q 


2 W f B\ 2 
Vb) 


n ( nr 0 + nr + r) (r — r 0 ) 1+1/n 
(2 n 2 + 3n + 1) 


tb 


(115) 


(116) 


Q- ( B X 

n (nr Q + nr B + r B ) (r B - r 0 ) 1+1/n 


(2n 2 + 3?r + 1) j 


(117) 


When n — 1, with C = C . the formula reduces to the Bingham, Equation 113; 
when t 0 = 0, with C = k, the formula reduces to the power law, Equation 84; 
and when n — 1 and r G = 0, with C = /x 0 , the formula reduces to the Newtonian, 
Equation 80, as it should be. 
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